#include <algorithm/analytical_model/algorithm/cpu/ConvolutionMethod.h>
#include <algorithm/analytical_model/algorithm/cpu/VectorConvert.h>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/Data.h>

using namespace solar::cpu;
using namespace solar::cuda::inter_op;

auto solar::cpu::intCauchyPeak(const Float3& x_axis, const Float3& y_axis,
                               const Float3* const interv, const float& s, const float& area,
                               float* gl_params, const Float3& focus_center, const float& sigma1)
    -> float
{
    const int m = 6;
    const int n = 6;
    const int gl_p1 = m * (m - 1);
    const int gl_p2 = n * (n - 1);
    float sum = 0;
    Vec3 interv_vec[3];
    auto x_axis_vec = toVec3(x_axis);
    auto y_axis_vec = toVec3(y_axis);
    auto focus_center_vec = toVec3(focus_center);
    for (int i = 0; i < 3; i++)
    {
        interv_vec[i] = toVec3(interv[i]);
    }
    for (int i = 0; i < m; i++)
    {
        float tmp_sum = 0;
        const float xi = gl_params[gl_p1 + i];
        const float ai = gl_params[gl_p1 + m + i];

        auto start_pt = focus_center_vec + (interv_vec[1] - interv_vec[0]) / 2 * xi;
        for (int j = 0; j < n; j++)
        {
            float xj = gl_params[gl_p2 + j];
            float aj = gl_params[gl_p2 + n + j];
            auto u = solar::dot(start_pt - focus_center_vec, x_axis_vec) +
                     xj * (interv_vec[2] - interv_vec[1]).length() / 2;
            auto v = solar::dot(start_pt - focus_center_vec, y_axis_vec);
            const auto tmp_num = 1 + sigma1 * sigma1 * (u * u + v * v);
            tmp_sum += aj * sigma1 * sigma1 / kMathPi / (tmp_num * tmp_num);
        }
        sum += tmp_sum * ai;
    }
    return sum * area / 4;
}
